MODULE GLOBAL
USE OMP_LIB
    
IMPLICIT NONE 

REAL(8), PARAMETER::mu=1d0,alpha=0.75d0,r=0.02d0,sigma=2d0,sigmag=2d0,omegac=0.3d0  
REAL(8), PARAMETER::pm=1d0,gamma=0.43d0,yloss=0.02d0,phi=0.18d0
!baseline,high/low beta, low kappa
REAL(8), PARAMETER::rhoa=0.777d0,sigmaa=0.029d0,mua=0d0,zt=0d0,taxrate=0.19d0
!lumpsum taxes
!REAL(8), PARAMETER::rhoa=0.777d0,sigmaa=0.029d0,mua=0d0,zt=0d0,taxrate=2d0
!Low tax
!REAL(8), PARAMETER::rhoa=0.777d0,sigmaa=0.029d0,mua=0d0,zt=0d0,taxrate=0.17d0
!high tax:
!REAL(8), PARAMETER::rhoa=0.777d0,sigmaa=0.029d0,mua=0d0,zt=0d0,taxrate=0.21d0
REAL(8), PARAMETER::gnmin=0d0,gnmax=1d0-1d-6,hmin=1d-6,hmax=1d0,width=1.5d0,in=1d0
REAL(8), PARAMETER::lambda=0.184d0,coup=0.0d0,pi_number=3.1415926535897932d0
!10-yr maturity
!REAL(8), PARAMETER::lambda=0.082d0,coup=0.0d0,pi_number=3.1415926535897932d0
REAL(8), PARAMETER::kappaval=1000000d0

!Rep agent:
!REAL(8),PARAMETER::alphac=1d0
!Low kappa:
!REAL(8),PARAMETER::alphac=0.4d0
!all else
REAL(8),PARAMETER::alphac=0.7d0
REAL(8)::constant_gn,transf_global
INTEGER,PARAMETER::welfare_case=0,only_repay=1
    
REAL(8), PARAMETER::tol=1d-6,weightq=1d0
!low/high beta
!INTEGER(4), PARAMETER::ns=21,nk=1,na=ns*nk,nb=71,nn=101,nmc=10000,n0=1000,nb2=501,na2=21,ndevs=3
!!all else
INTEGER(4), PARAMETER::ns=21,nk=1,na=ns*nk,nb=61,nn=101,nmc=10000,n0=1000,nb2=501,na2=21,ndevs=3
INTEGER(4), PARAMETER::nt=9000,nt2=32,nsam=1000,nm=100,nf=2    
INTEGER(4), PARAMETER::maxiter=4000,niterv=15

!baseline
INTEGER(4), PARAMETER::ncdf=70,nquad_v=15,neps= nquad_v,nquad_q=11
!high beta
!INTEGER(4), PARAMETER::ncdf=70,nquad_v=11,neps= nquad_v,nquad_q=11

INTEGER(4), PARAMETER::save_results=1,choice_defcost=2,choice_counter=0,choice_coup=0,choice_bound=1
INTEGER(4), PARAMETER::choice_dist=0,choice_nog=0,choice_fe =0,choice_nowbar=0      !for choice_nowbar=1 set wbar=0d0
INTEGER(4), PARAMETER::choice_correl=0,choice_tax_hat=0                             !=0 no tax_hat, =1 frictionless tax level
INTEGER(4), PARAMETER::choice_rule=0                                                !=0 use saved data, =1 use given coefficient for gN/yN
INTEGER(4),PARAMETER::nnaux=choice_fe*nn+(1-choice_fe)*1,export_data=0              !=1 to export subsamples series
INTEGER(4),PARAMETER::choice_constant_gn=0,choice_ols=0

REAL(8),PARAMETER::th_1=-0.27982d0,th_2=-1.21d0,mingn=0.045844d0,scale=1.01d0
REAL(8), PARAMETER::cdfmin=3.167d-5,cdfmax=1d+0-3.167d-5,huge_n=-1d6
REAL(8), PARAMETER::resp_coefl1=0.16d0,resp_coefl2=0d0*0.4248d0,resp_coefl3=0d0*0.539928d0,resp_coefl4=0d0*(lambda+r)
REAL(8), PARAMETER::resp_coefh1=0.16d0,resp_coefh2=-0d0*0.297183d0,resp_coefh3=-0d0*0.257716d0,resp_coefh4=0d0*(lambda+r)
REAL(8), PARAMETER::resp_coef1=0.200627d0,resp_coef2=-0.115658d0,resp_coef3=-0.026376d0*(lambda+r),resp_coef4=0d0*(lambda+r)
REAL(8)::beta,wbar,icost,chi,psi1,psi2,rate_time
REAL(8)::agrid(na),bgrid(nb),bgrid2(nb2),gngrid(nn),hgrid(nn),defaultgrid(2),prob(na2,na2),m2(10),amin,amax,aautgrid(na),udefgrid(na)
REAL(8)::kw,aux1,aux3,aux4,aux5,meanb2
REAL(8)::v0(na2,nb2),v1(na2,nb2),vaut0(na2),vaut1(na2),vcon0(na2,nb2),vcon1(na2,nb2),q1(na2,nb2)
REAL(8), DIMENSION(na2,nb2)::p1,h1,gn1,tau1,ct1,w1
REAL(8), DIMENSION(na2)::paut1,tauaut1,haut1,gnaut1,ctaut1,waut1

REAL(8)::cdfgrid(ncdf),epsgrid(neps),epsmin,epsmax,quad_w_v(1:nquad_v), quad_x_v(1:nquad_v), quad_w_hermite(1:neps),&
quad_x_hermite(1:neps), quad_w_q(1:nquad_q), quad_x_q(1:nquad_q)
REAL(8)::break_matrix(na,nb), break_matrix_q(na,nb)
REAL(8)::coeff_matrix(na,4,nb), coeff_matrix_q(na,4,nb),v0_matrix(nb,na), vaut_matrix(nb,na) 
REAL(8)::BREAK_eps(ncdf),CSCOEF_eps(4,ncdf)
REAL(8)::break_matrix_w(na,nb), coeff_matrix_w(na,4,nb),bnext1(na2,nb2)

REAL(8)::agrid2(na2)

REAL(8), DIMENSION(nb,na)::v_matrix, w_matrix,default_decision,b_next_matrix,x_matrix_new,gn_matrix_new,v0_matrixA_w
REAL(8), DIMENSION(nb,na)::b0_next_matrix, b1_next_matrix, q_matrix, q_matrix_nodef
REAL(8)::wgloboff(nn),wautgloboff(nn),gnnglob(nn),payoff

REAL(8)::bmin,bmax,reg_global,x_global,x_global1,x_global2 ,x_global3,gn_global,gn_global1,gn_global2,gn_global3,tax_hat,ctrad_global,ct_global1,ct_global2,ct_global3,bp_global
  
INTEGER::bp1(na2,nb2),nout,choice_glob,counter_global,monot_v0,index_bfeas2(na),def1(na2,nb2)
INTEGER::start_time    
INTEGER::seed,hours,mins,ioff_globinfeasdebt(na),iterfun,i_bzero,i_bzero2,version_code,index_bfeas4(na)

REAL(8)::h_base(na,nb), wage_base(na,nb), gn_base(na,nb),transf_base(na,nb),bp_base(na,nb), haut_base(na), gnaut_base(na),  v_base(nb, na),v0_base(nb, na), vaut_base(nb, na)   
REAL(8)::q1base(na,501),bgridbase(501)

REAL(8),DIMENSION(nb,na)::b_next_matrix_2,q_matrix_2,q_matrix_20,b_next_matrix_20
REAL(8),DIMENSION(nb2,na2)::gn1_2,h1_2,ct1_2,q1_2,q1_20,gn1_20,h1_20,ct1_20,vcon1_2,vcon1_20,v1_2,v1_20,b1_2,b1_20
REAL(8),DIMENSION(na2)::gnaut1_2,ctaut1_2,haut1_2,gnaut1_20,ctaut1_20,haut1_20,vaut1_2,vaut1_20
INTEGER,DIMENSION(nb,na)::default_decision_2,default_decision_20
INTEGER,DIMENSION(nb2,na2)::def1_2,def1_20,bp1_2,bp1_20
INTEGER::index_bfeasmax2(na)
REAL(8)::q_nodef1(na2,nb2)

!!!!$OMP NOUSETHREADPRIVATE(a_initial,a0,aux0,b_initial,b0,ct_global,ftype,x_global,x_global1,x_global2,reg_global,gn_global,gn_global1,gn_global2,b_global,i_default_global,choice_glob,rate_time,start_time)
    
CONTAINS

!Implement Tauchen (1986) Method
SUBROUTINE BUILD_TRANSITION(N,h,rho_h,sigma_h,alpha_h,pi)

IMPLICIT NONE

INTEGER(4),INTENT(IN)::N
REAL(8),INTENT(IN)::h(N),rho_h,sigma_h,alpha_h
REAL(8),INTENT(OUT)::pi(N,N)
REAL(8)::h_sort(N),temp1,w,temp1_low,temp1_high
INTEGER(4) I_h(N),i,count,k

h_sort = h

w =  h_sort(2) - h_sort(1)

DO i=1,N
	temp1 = (h_sort(1) - rho_h*h_sort(i) - alpha_h + 0.5d0*w)/sigma_h
	pi(i,1) = alnorm(temp1)
ENDDO

w = h_sort(N) - h_sort(N-1)

DO i=1,N
	temp1 = (h_sort(N) - rho_h*h_sort(i) - alpha_h - 0.5d0*w)/sigma_h
	pi(i,N) = 1d0 - alnorm(temp1)	
ENDDO

DO k=2,N-1
	w = h_sort(k) - h_sort(k-1)

	DO i=1,N
	
	temp1_high = (h_sort(k) - rho_h*h_sort(i) - alpha_h + 0.5d0*w)/sigma_h
	temp1_low = (h_sort(k) - rho_h*h_sort(i) - alpha_h - 0.5d0*w)/sigma_h
	pi(i,k) = alnorm(temp1_high) - alnorm(temp1_low)
	ENDDO
ENDDO

END SUBROUTINE BUILD_TRANSITION
!*********************************************************************************
FUNCTION LOCATE(xx,x)

IMPLICIT NONE
      
REAL(8),INTENT(IN)::xx(:)
REAL(8),INTENT(IN)::x
INTEGER::locate
INTEGER::n,jl,jm,ju
LOGICAL::ascnd
      
n=size(xx)
ascnd = (xx(n) >= xx(1))
      
jl=0
ju=n+1

DO
    IF (ju-jl <= 1) EXIT
         
    jm=(ju+jl)/2
    IF (ascnd .eqv. (x >= xx(jm))) THEN
        jl=jm
    ELSE
        ju=jm
    ENDIF
ENDDO
      
      
IF (x <= xx(1)) THEN
    locate=1
ELSEIF (x >= xx(n)) THEN
    locate=n
ELSE
    locate=jl
ENDIF
      
END FUNCTION LOCATE
!*********************************************************************************
FUNCTION LOCATE2(xx,x)

IMPLICIT NONE
      
REAL(8),INTENT(IN)::xx(:)
REAL(8),INTENT(IN)::x
INTEGER::locate2
INTEGER::n,jl
      
n=size(xx)
jl=locate(xx,x)
locate2=jl

IF (DABS(xx(jl)-x).GT.DABS(xx(MIN(n,jl+1))-x)) THEN
    locate2=MIN(n,jl+1)
ELSEIF (DABS(xx(jl)-x).GT.DABS(xx(MAX(1,jl-1))-x)) THEN
    locate2=MAX(1,jl-1)
ENDIF

END FUNCTION LOCATE2
!*********************************************************************************
REAL FUNCTION LINEAR_INT(x,y,z)

IMPLICIT NONE
 
REAL(8),INTENT(IN)::x(:),y(:),z
INTEGER::i,n
      
n=SIZE(x)
      
IF (n/=SIZE(y)) THEN
    WRITE(0,*)'ERROR IN LINEAR INTERPOLATION: X AND Y MUST AGREE'
    STOP
ENDIF
      
IF (z<x(1)) THEN
    i=1
ELSEIF (z>=x(n)) THEN
    i=n-1
ELSE
    i=locate(x,z)
ENDIF
      
linear_int=y(i)+(y(i+1)-y(i))/(x(i+1)-x(i))*(z-x(i))
      
END FUNCTION LINEAR_INT
!*********************************************************************************
SUBROUTINE INITIAL

IMPLICIT NONE

INTEGER::i,j,i_a,i_b,ii_b,nbaux,index_opt,ILEFT, IRIGHT,NINTV,i0(1)
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb)
REAL(8)::new_value
EXTERNAL DCSDEC

v0_matrix     = 0d0
vaut_matrix   = -1d4
q_matrix_nodef= payoff/ (lambda + r)

index_bfeas4 = 0d0
    
!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
ENDDO
index_bfeas2=MAXVAL(index_bfeas4)

do i_a = 1,na                 
    IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
        ii_b=index_bfeas4(i_a)+1
        nbaux =nb-ii_b+1
        FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
	    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	    DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	    call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
    ENDIF         
ENDDO
         
do i_a = 1,na
    FDATA = q_matrix_nodef(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
    break_matrix_q(i_a, :) = BREAK_GRID
    coeff_matrix_q(i_a, :,:) = CSCOEF
end do


PRINT*,'Initial Values for solving model uploaded'
PRINT*,' '

END SUBROUTINE
!*********************************************************************************
SUBROUTINE INITIAL_POL

IMPLICIT NONE

INTEGER::i,j,i_a,i_b,ii_b,nbaux,index_opt,ILEFT, IRIGHT,NINTV,i0(1)
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb)
REAL(8)::new_value
EXTERNAL DCSDEC

!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt
ENDDO

index_bfeas2=MAXVAL(index_bfeas4)

do i_a = 1,na         
    IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
        ii_b=index_bfeas4(i_a)+1
        nbaux =nb-ii_b+1
        FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
	    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	    DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	    call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
    ENDIF
         
ENDDO

do i_a = 1,na
    FDATA = q_matrix_nodef(:,i_a)
    ILEFT  = 0
    IRIGHT = 0
    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
    DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
    call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
    break_matrix_q(i_a, :) = BREAK_GRID
    coeff_matrix_q(i_a, :,:) = CSCOEF
end do
     
OPEN(1,FILE='output//gvcon.txt')
OPEN(2,FILE='output//gp.txt')
OPEN(3,FILE='output//gh.txt')
OPEN(4,FILE='output//gibp.txt')
OPEN(5,FILE='output//gct.txt')
OPEN(7,FILE='output//ggn.txt')
OPEN(8,FILE='output//gtau.txt')    
OPEN(9,FILE='output//gw.txt')    
OPEN(10,FILE='output//gq.txt')    
OPEN(11,FILE='output//gdef1.txt')    
            
DO i = 1,na2
DO j = 1,nb2
    READ(1,*) vcon1(i,j)
    READ(2,*) p1(i,j)
    READ(3,*) h1(i,j)
    READ(4,'(I7)') bp1(i,j)
    READ(5,*) ct1(i,j)
    READ(7,*) gn1(i,j)
    READ(8,*) tau1(i,j)
    READ(9,*) w1(i,j)
    READ(10,*) q1(i,j)
    READ(11,*) def1(i,j)
ENDDO
ENDDO
    
CLOSE(1)
CLOSE(2)
CLOSE(3)
CLOSE(4)
CLOSE(5)
CLOSE(7)
CLOSE(8)
CLOSE(9)
CLOSE(10)
CLOSE(11)
    
OPEN(1,FILE='output//gvaut.txt')
OPEN(2,FILE='output//gpaut.txt')
OPEN(3,FILE='output//ghaut.txt')
OPEN(5,FILE='output//gctaut.txt')
OPEN(7,FILE='output//ggnaut.txt')
OPEN(8,FILE='output//gtauaut.txt')    
OPEN(9,FILE='output//gwaut.txt')    
        
DO i = 1,na2
    READ(1,*) vaut1(i)
    READ(2,*) paut1(i)
    READ(3,*) haut1(i)
    READ(5,*) ctaut1(i)
    READ(7,*) gnaut1(i)
    READ(8,*) tauaut1(i)
    READ(9,*) waut1(i)
ENDDO
    
CLOSE(1)
CLOSE(2)
CLOSE(3)
CLOSE(5)
CLOSE(7)
CLOSE(8)
CLOSE(9)

OPEN(1,FILE='output//bgrid2.txt')    
            

DO i=1,nb2
    READ(1,*) bgrid2(i)
    v1(:,i)=MAX(vcon1(:,i),vaut1(:))
    def1(:,i)=-(vcon1(:,i).GE.vaut1(:))
ENDDO

CLOSE(1)

i0 = LOCATE2(bgrid2,0d0)
i_bzero2 = i0(1)
bgrid2(i_bzero2)=0d0

OPEN(1,FILE='output//agrid2.txt')           
DO i = 1,na2
    READ(1,*) agrid2(i)
ENDDO    
CLOSE(1)

CALL BUILD_TRANSITION(na2,agrid2,rhoa,sigmaa,mua,prob)

PRINT*,'Solution for policy fcns and value fcns uploaded'
PRINT*,' '

END SUBROUTINE
!*********************************************************************************
SUBROUTINE INITIAL_V

IMPLICIT NONE

INTEGER::i,j,i_a,i_b,ii_b
   
OPEN(1,FILE='output//gvcon.txt')
OPEN(2,FILE='output//gq_nodef.txt')
            
DO i = 1,na2
DO j = 1,nb2
    READ(1,*) vcon1(i,j)
    v0_matrix(j,i) = vcon1(i,j)
    READ(2,*) q_matrix_nodef(j,i) 
ENDDO
ENDDO
    
CLOSE(1)
CLOSE(2)
    
OPEN(1,FILE='output//gvaut.txt')
        
DO i = 1,na2
    READ(1,*) vaut1(i)
    vaut_matrix(:,i)=vaut1(i)
ENDDO
    
CLOSE(1)


PRINT*,'Solution for value fcns and prices uploaded'
PRINT*,' '

END SUBROUTINE INITIAL_V
!*********************************************************************************
SUBROUTINE TIME_SERIES(Ntime,Nstates,Pi,ystate,y0,ind_state1)

IMPLICIT NONE

INTEGER(4),INTENT(IN)::Ntime,Nstates
REAL(8),INTENT(IN)::Pi(Nstates,Nstates),y0,ystate(Nstates)   
INTEGER(4),INTENT(OUT)::ind_state1(Ntime)
REAL(8)::dummy11(1,1),dum3(Nstates),ind_state_prime
INTEGER(4)::i,j,ind_state,dummy2(1)


dummy2 = MINLOC( (y0 - ystate)**2d0)

ind_state = INT(dummy2(1))
ind_state1(1) = INT(dummy2(1))

!Start w y0 or closest to y0

DO i=2,Ntime
	dum3(1) = Pi(ind_state,1)

	DO j=2,Nstates
		dum3(j) = dum3(j-1) + Pi(ind_state,j)
	ENDDO

	CALL r8vec_uniform_01(1,seed,dummy11)
    
    IF (dum3(Nstates).GT.1d0-1d-5) dum3(Nstates)=1d0

	ind_state1(i) = 0

	IF ( dummy11(1,1) .LE. dum3(1) ) THEN
		ind_state1(i) = 1
	ELSE
		DO j=1,Nstates-1
		IF (dum3(j+1) .GT. dum3(j)) THEN

			IF ( (dummy11(1,1) .GT. dum3(j)) .AND. (dummy11(1,1) .LE. dum3(j+1)) ) THEN
				ind_state1(i) = j+1
			ENDIF
		ENDIF
		ENDDO
    ENDIF
   
	ind_state = ind_state1(i)
    
    IF (ind_state1(i).EQ.0) THEN
        PRINT*,'Problem w/TIME_SERIES',dum3
        PRINT*,dummy11
    PAUSE
    STOP
    ENDIF
ENDDO

END SUBROUTINE TIME_SERIES
!*********************************************************************************
FUNCTION mean(x)

IMPLICIT NONE

REAL(8),DIMENSION(:),INTENT(IN)::x
REAL(8)::mean
INTEGER::n

n = SIZE(x)

mean = SUM(x)/DBLE(n)

END FUNCTION mean
!*********************************************************************************
FUNCTION std(x)

IMPLICIT NONE

REAL(8),DIMENSION(:),INTENT(IN)::x
REAL(8)::std,mm
INTEGER::n

n = SIZE(x)
mm = mean(x)

std = DSQRT(SUM((x-mm)**2d0)/DBLE(n))

END FUNCTION std
!*********************************************************************************
FUNCTION corr(x,xx)

IMPLICIT NONE

REAL(8),DIMENSION(:),INTENT(IN)::x,xx
REAL(8)::corr,aux,m1,m2,s1,s2
INTEGER::n

n = SIZE(x)
m1 = mean(x)
m2 = mean(xx)
s1 = std(x)
s2 = std(xx)

IF (s1*s2 .EQ. 0d0) THEN
!    PRINT*, 'Problem Correlation: stddev<=0',s1,s2
!    STOP
ENDIF

aux = (SUM((x-m1)*(xx-m2))/DBLE(n))

IF (s1*s2 .GT. 0d0) THEN
    corr = aux/(s1*s2)
ELSE
    corr=0d0
ENDIF

END FUNCTION corr
!**************************************************************
!* This function computes the quantiles for a time series     *
!**************************************************************
FUNCTION quantile(x,p)
! Calculate Quartiles using SAS Method 5
! This method is the default method of SAS and is based on the empirical distribution function.
! Based on discussion in this paper http://www.haiweb.org/medicineprices/manual/quartiles_iTSS.pdf

IMPLICIT NONE
REAL(8), DIMENSION(:),INTENT(IN)::x(:)
REAL(8), ALLOCATABLE::x2(:),xx(:)
REAL(8), INTENT(IN)::p
REAL(8)::temp,ii,d,qq,quantile
INTEGER::n,i,iu,il,j

   n = SIZE(x);

   ALLOCATE(xx(n+2),x2(n))

   x2=x;
   
   DO j=1,n
   DO i=1,n-1
    IF (x2(i) > x2(i+1)) THEN
        temp=x2(i)
        x2(i)=x2(i+1)
        x2(i+1)=temp
    ENDIF
   ENDDO
   ENDDO
      
   xx = [x2(1), x2, x2(n)];  
   
   ii = p*DBLE(n)+1.5d0;
   iu = CEILING(ii);
   il = FLOOR(ii);
   d = (ii-il);
   quantile = xx(il)*(1d0-d)+xx(iu)*d;

END FUNCTION quantile
!******************************************************************************************
!* This routine generates a random number from a Uniform [0,1] distribution
!******************************************************************************************
subroutine r8vec_uniform_01 ( n, seed1, r )

!! R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    31 May 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Second Edition,
!    Springer, 1987,
!    ISBN: 0387964673,
!    LC: QA76.9.C65.B73.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, December 1986, pages 362-376.
!
!    Pierre L'Ecuyer,
!    Random Number Generation,
!    in Handbook of Simulation,
!    edited by Jerry Banks,
!    Wiley, 1998,
!    ISBN: 0471134031,
!    LC: T57.62.H37.
!
!    Peter Lewis, Allen Goodman, James Miller,
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, 1969, pages 136-143.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of entries in the vector.
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, which
!    should NOT be 0.  On output, SEED has been updated.
!
!    Output, real ( kind = 8 ) R(N), the vector of pseudorandom values.
!

  implicit none


  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ) k
  integer ( kind = 4 ) seed1
  real    ( kind = 8 ) r(n)

  seed = seed1
  
  do i = 1, n

    k = seed / 127773

    seed = 16807 * ( seed - k * 127773 ) - k * 2836

    if ( seed < 0 ) then
      seed = seed + 2147483647
    end if

    r(i) = real ( seed, kind = 8 ) * 4.656612875D-10

  end do

  return

end subroutine r8vec_uniform_01
!*********************************************************************************
FUNCTION alnorm(x)

!! ALNORM computes the cumulative density of the standard normal distribution.
!
!  Modified:
!
!    13 January 2008
!
!  Author:
!
!    Original FORTRAN77 version by David Hill.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    David Hill,
!    Algorithm AS 66:
!    The Normal Integral,
!    Applied Statistics,
!    Volume 22, Number 3, 1973, pages 424-427.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, is one endpoint of the semi-infinite interval
!    over which the integration takes place.
!
!    Input, logical UPPER, determines whether the upper or lower
!    interval is to be integrated:
!    .TRUE.  => integrate from X to + Infinity;
!    .FALSE. => integrate from - Infinity to X.
!
!    Output, real ( kind = 8 ) ALNORM, the integral of the standard normal
!    distribution over the desired interval.
!
  implicit none

  real ( kind = 8 ), parameter :: a1 = 5.75885480458D+00
  real ( kind = 8 ), parameter :: a2 = 2.62433121679D+00
  real ( kind = 8 ), parameter :: a3 = 5.92885724438D+00
  real ( kind = 8 ) alnorm
  real ( kind = 8 ), parameter :: b1 = -29.8213557807D+00
  real ( kind = 8 ), parameter :: b2 = 48.6959930692D+00
  real ( kind = 8 ), parameter :: c1 = -0.000000038052D+00
  real ( kind = 8 ), parameter :: c2 = 0.000398064794D+00
  real ( kind = 8 ), parameter :: c3 = -0.151679116635D+00
  real ( kind = 8 ), parameter :: c4 = 4.8385912808D+00
  real ( kind = 8 ), parameter :: c5 = 0.742380924027D+00
  real ( kind = 8 ), parameter :: c6 = 3.99019417011D+00
  real ( kind = 8 ), parameter :: con = 1.28D+00
  real ( kind = 8 ), parameter :: d1 = 1.00000615302D+00
  real ( kind = 8 ), parameter :: d2 = 1.98615381364D+00
  real ( kind = 8 ), parameter :: d3 = 5.29330324926D+00
  real ( kind = 8 ), parameter :: d4 = -15.1508972451D+00
  real ( kind = 8 ), parameter :: d5 = 30.789933034D+00
  real ( kind = 8 ), parameter :: ltone = 7.0D+00
  real ( kind = 8 ), parameter :: p = 0.398942280444D+00
  real ( kind = 8 ), parameter :: q = 0.39990348504D+00
  real ( kind = 8 ), parameter :: r = 0.398942280385D+00
  logical up
  logical upper
  real ( kind = 8 ), parameter :: utzero = 18.66D+00
  real ( kind = 8 ) x
  real ( kind = 8 ) y
  real ( kind = 8 ) z

  
  upper = .FALSE.

  up = upper
  z = x

  if ( z < 0.0D+00 ) then
    up = .not. up
    z = - z
  end if

  if ( ltone < z .and. ( ( .not. up ) .or. utzero < z ) ) then

    if ( up ) then
      alnorm = 0.0D+00
    else
      alnorm = 1.0D+00
    end if

    return

  end if

  y = 0.5D+00 * z * z

  if ( z <= con ) then

    alnorm = 0.5D+00 - z * ( p - q * y &
      / ( y + a1 + b1 &
      / ( y + a2 + b2 &
      / ( y + a3 ))))

  else

    alnorm = r * exp ( - y ) &
      / ( z + c1 + d1 &
      / ( z + c2 + d2 &
      / ( z + c3 + d3 &
      / ( z + c4 + d4 &
      / ( z + c5 + d5 &
      / ( z + c6 ))))))

  end if

  if ( .not. up ) then
    alnorm = 1.0D+00 - alnorm
  end if

  return

end function alnorm
!**************************************************************************
! Test function for Brent method
REAL(8) FUNCTION AFunction(ftype2,a02,aux02,ct_global2,x)

IMPLICIT NONE

REAL(8),INTENT(IN)::a02,aux02,ct_global2
INTEGER,INTENT(IN)::ftype2
REAL(8)::x,aahat,aa1,dhcn,ucn,ucnaux1,ucnaux2,uaux
REAL(8)::lambdatemp,temp1,temp2,ptemp,gnaux,caux,ugnaux,c1aux,paux,taux,haux,cnaux,ctaux

IF (ftype2.EQ. 35) THEN
    haux = (wbar/alpha/x)**(1d0/(alpha-1))
    gnaux = aux02
    ctaux = a02+x*gnaux-taxrate*(a02+in*x*haux**alpha)
    cnaux = haux**alpha-gnaux
    AFunction = x-(1d0-omegac)/omegac*(ctaux/cnaux)**(1d0+mu)   
ELSEIF (ftype2.EQ. 36) THEN
    ctaux = ct_global2
    gnaux = aux02
    paux  = (1d0-omegac)/omegac*(ctaux/(x**alpha-gnaux))**(1d0+mu)
    AFunction = alpha*paux*(x**(alpha-1d0))-wbar
ELSEIF (ftype2 .EQ. 100) THEN
    paux = (1d0-omegac)/omegac*(ct_global2/x)**(1d0+mu)
    IF (mu.EQ.0d0) THEN
        caux = (ct_global2**omegac)*(x**(1d0-omegac))
        c1aux = (1d0-omegac)*(cT_global2**omegac)*(x**(-omegac))    
    ELSE
        caux = (omegac*ct_global2**(-mu)+(1d0-omegac)*x**(-mu))**(-1d0/mu)
        c1aux = (1d0-omegac)*((omegac*ct_global2**(-mu)+(1d0-omegac)*x**(-mu))**(-1d0/mu-1d0))*(x**(-mu-1d0))
    ENDIF
    AFunction = chi*(1d0-x)**(-sigmag) - (1d0-chi)*(caux**(-sigma))*c1aux
ELSEIF (ftype2 .EQ. 101) THEN
        paux = (1d0-omegac)/omegac*(ct_global2/x)**(1d0+mu)
    temp1 = (wbar/alpha/paux)**(1d0/(alpha-1d0))  
    temp1 = (alpha/wbar*paux)**(1d0/(1d0-alpha))  
    gnaux = temp1**alpha-x  
    ugnaux = chi*(gnaux**(-sigmag))
    IF (mu.EQ.0d0) THEN
        caux = (ct_global2**omegac)*(x**(1d0-omegac))
        c1aux = (1d0-omegac)*(cT_global2**omegac)*(x**(-omegac))    
    ELSE
        caux = (omegac*ct_global2**(-mu)+(1d0-omegac)*x**(-mu))**(-1d0/mu)
        c1aux = (1d0-omegac)*((omegac*ct_global2**(-mu)+(1d0-omegac)*x**(-mu))**(-1d0/mu-1d0))*(x**(-mu-1d0))
    ENDIF
    temp1=MAX(temp1,-1d6)
    temp1=MIN(temp1,1d6)
    aahat  = 1d0/(temp1+(1d0-temp1)*alphac)
    aa1    = (temp1+(1d0-temp1)*alphac**(1d0-sigma))
    IF (alphac.EQ.1d0) THEN 
        aa1=1d0
        aahat = 1d0
    ENDIF    
    uaux=(caux**(1d0-sigma))/(1d0-sigma)
    ucn = caux**(-sigma)*c1aux*(aahat**(1d0-sigma))*aa1
    dhcn= (1d0+mu)*temp1/(alpha-1d0)/x
    ucnaux1 = (aahat**(1d0-sigma))*(1d0-alphac**(1d0-sigma))*dhcn*uaux
    ucnaux2 = (1d0-sigma)*(aahat**(-sigma))*(-1)*(1d0-alphac)/((temp1+(1d0-temp1)*alphac)**2d0)*dhcn*aa1*uaux
    AFunction = (1d0-chi)*(ucn + ucnaux1 +ucnaux2) + ugnaux*(alpha*temp1**(alpha-1d0)*dhcn-1d0)
    RETURN
ELSEIF (ftype2 .EQ. 102) THEN
    paux = (1d0-omegac)/omegac*(ct_global2/x)**(1d0+mu)
    taux = paux*(1d0-x)-aux02
    IF (mu.EQ.0d0) THEN
        caux = (ct_global2**omegac)*(x**(1d0-omegac))
        c1aux = (1d0-omegac)*(cT_global2**omegac)*(x**(-omegac))    
    ELSE
        caux = (omegac*ct_global2**(-mu)+(1d0-omegac)*x**(-mu))**(-1d0/mu)
        c1aux = (1d0-omegac)*((omegac*ct_global2**(-mu)+(1d0-omegac)*x**(-mu))**(-1d0/mu-1d0))*(x**(-mu-1d0))
    ENDIF
    AFunction = chi*(1d0-x)**(-sigmag) - (1d0-chi)*(caux**(-sigma))*c1aux-icost*2d0*(taux-tax_hat)*paux*((1d0+mu)*(1d0-x)/x+1d0)
ELSEIF (ftype2 .EQ. 103) THEN
    paux = (1d0-omegac)/omegac*(ct_global2/x)**(1d0+mu)
    temp1 = (wbar/alpha/paux)**(1d0/(alpha-1d0))  !ok
    taux = paux*(temp1**alpha-x)-aux02
    gnaux = temp1**alpha-x  !ok
    ugnaux = chi*(gnaux**(-sigmag))
    IF (mu.EQ.0d0) THEN
        caux = (ct_global2**omegac)*(x**(1d0-omegac))
        c1aux = (1d0-omegac)*(cT_global2**omegac)*(x**(-omegac))    
    ELSE
        caux = (omegac*ct_global2**(-mu)+(1d0-omegac)*x**(-mu))**(-1d0/mu)
        c1aux = (1d0-omegac)*((omegac*ct_global2**(-mu)+(1d0-omegac)*x**(-mu))**(-1d0/mu-1d0))*(x**(-mu-1d0))
    ENDIF
    AFunction = (1d0-chi)*(caux**(-sigma))*c1aux+ ugnaux*((temp1**alpha)*alpha*(1d0+mu)/(alpha-1d0)/x-1d0)+icost*2d0*(taux-tax_hat)*paux*((temp1**alpha)*(1d0+mu)/x/(1d0-alpha)-mu)
ELSEIF (ftype2 .EQ. 198) THEN
    gnaux = aux02 
    paux  = (1d0-omegac)/omegac*(ct_global2/((1d0-gnaux)*x**alpha))**(1d0+mu)
    AFunction = alpha*paux*(x**(alpha-1d0))-wbar

ELSEIF (ftype2 .EQ. 199) THEN
    gnaux = aux02 
    paux  = (1d0-omegac)/omegac*(ct_global2/(x**alpha-gnaux))**(1d0+mu)
    AFunction = alpha*paux*(x**(alpha-1d0))-wbar

ELSEIF (ftype2.EQ.200) THEN
    gnaux     = aux02
    AFunction = (1d0-gnaux)*(kw*x)**(1d0/(1d0+mu))-x*gnaux - (a02-taxrate*(a02+in*x))
ELSEIF (ftype2.EQ.302) THEN
    paux  = (1d0-omegac)/omegac*(ct_global2/x)**(1d0+mu)
    haux  = (alpha/wbar*paux)**(1d0/(1d0-alpha))
    AFunction = paux*((1d0-in*taxrate)*haux**alpha-x) - taxrate*a02 -aux02 
ELSEIF (ftype2 .EQ. 401) THEN
    cnaux = aux02
    paux = (1d0-omegac)/omegac*(x/cnaux)**(1d0+mu)
    AFunction = x-a02+taxrate*(a02+paux*in)-paux*(1d0-cnaux)
ELSEIF (ftype2 .EQ. 402) THEN
    paux = (1d0-omegac)/omegac*(ct_global2/x)**(1d0+mu)
    AFunction = paux*(1d0-x-in*taxrate)-aux02-taxrate*a02
ELSEIF (ftype2 .EQ. 405) THEN
    paux = (1d0-omegac)/omegac*(ct_global2/x)**(1d0+mu)
    IF (mu.EQ.0d0) THEN
        caux = (ct_global2**omegac)*(x**(1d0-omegac))
        c1aux = (1d0-omegac)*(cT_global2**omegac)*(x**(-omegac))    
    ELSE
        caux = (omegac*ct_global2**(-mu)+(1d0-omegac)*x**(-mu))**(-1d0/mu)
        c1aux = (1d0-omegac)*((omegac*ct_global2**(-mu)+(1d0-omegac)*x**(-mu))**(-1d0/mu-1d0))*(x**(-mu-1d0))
    ENDIF
    ucn = (1d0-chi)*(caux**(-sigma))*c1aux
    gnaux  = ((aux02+taxrate*a02)/paux+taxrate*x)/(1d0-taxrate)
    haux   = (x+gnaux)**(1d0/alpha)
    ugnaux = chi*(gnaux**(-sigmag))
    lambdatemp = -ucn/(taxrate*paux+(gnaux-taxrate*haux**alpha)*(1d0+mu)*paux/x)
    AFunction = ugnaux-ucn-lambdatemp*(paux+(gnaux-taxrate*haux**alpha)*(1d0+mu)*paux/x)
ENDIF

END FUNCTION AFunction
!**********************************************************
! TRUE if x1*x2 negative
integer Function RootBracketed(x1,x2)
implicit none
  real(8) x1,x2 
  integer resultat
  if ((x1 > 0.and.x2 > 0).or.(x1 < 0.and.x2 < 0)) then 
    resultat = 0
  else
    resultat = 1
  endif
  RootBracketed = resultat
end Function RootBracketed

! returns the minimum of two real numbers
real(8) Function Minimum(x1,x2) 
  real(8) x1,x2,resultat
  if (x1 < x2) then
    resultat = x1
  else 
    resultat = x2
  endif
  Minimum = resultat
end Function Minimum
!*****************************************************
!*              Brent Method Function                *
!* ------------------------------------------------- *
!* The purpose is to find a real root of a real      *
!* function f(x) using Brent method.                 *
!*                                                   *
!* INPUTS:  x1,x2     : interval of root             *
!*          Tolerance : desired accuracy for root    *
!*          maxIter   : maximum number of iterations *
!*                                                   *
!* OUTPUTS: The function returns the root value      *
!*          ValueAtRoot : value of f(root)           *
!*          niter    : number of done iterations     *
!*          error    : =0, all OK                    *
!*                   : =1, no root found in interval *
!*                   : =2, no more iterations !      *
!*****************************************************
real(8) Function BrentRoots( ftype2,a02,aux02,ct_global2,x1, x2, Tolerance, maxIterations,&
		            valueAtRoot, niter, error)  

implicit none

REAL(8),INTENT(IN)::a02,aux02,ct_global2
INTEGER,INTENT(IN)::ftype2

  real(8),parameter::FPP = 1.d-11, nearzero = 1.d-20 

  real(8)  x1,x2,Tolerance,valueAtRoot  
  integer maxIterations,niter, error    

  real(8) resultat, AA, BB, CC, DD, EE, FA, FB, FC, Tol1, PP, QQ, RR, SS, xm
  integer i, done

  i = 0; done = 0;   error = 0
  AA = x1;  BB = x2;  FA = AFunction(ftype2,a02,aux02,ct_global2,AA); FB = AFunction(ftype2,a02,aux02,ct_global2,BB)
  if (RootBracketed(FA,FB).eq.0) then 
    error = 1
  else 
    FC = FB;
    do while (done.eq.0.and.i < maxIterations)
      if (RootBracketed(FC,FB).eq.0) then
        CC = AA; FC = FA; DD = BB - AA; EE = DD
      endif
      if (dabs(FC) < dabs(FB)) then
        AA = BB; BB = CC; CC = AA
        FA = FB; FB = FC; FC = FA
      endif
      Tol1 = 2.0 * FPP * dabs(BB) + 0.5 * Tolerance
      xm = 0.5 * (CC-BB)
      if ((dabs(xm) <= Tol1).or.(dabs(FA) < nearzero)) then
        ! A root has been found
        resultat = BB;
        done = 1
        valueAtRoot = AFunction(ftype2,a02,aux02,ct_global2,resultat)
      else 
        if ((dabs(EE) >= Tol1).and.(dabs(FA) > dabs(FB))) then
          SS = FB/ FA;
          if (dabs(AA - CC) < nearzero) then
            PP = 2.0 * xm * SS;
            QQ = 1.0 - SS;
          else 
            QQ = FA/FC;
            RR = FB /FC;
            PP = SS * (2.0 * xm * QQ * (QQ - RR) - (BB-AA) * (RR - 1.0));
            QQ = (QQ - 1.0) * (RR - 1.0) * (SS - 1.0);
          endif
          if (PP > nearzero) QQ = -QQ;
          PP = dabs(PP);
          if ((2.0 * PP) < Minimum(3.0*xm *QQ-dabs(Tol1 * QQ), dabs(EE * QQ))) then
            EE = DD;  DD = PP/QQ;
          else 
            DD = xm;   EE = DD;
          endif
        else 
          DD = xm;
          EE = DD;
        endif
        AA = BB;
        FA = FB;
        if (dabs(DD) > Tol1) then 
          BB = BB + DD;
        else 
          if (xm > 0) then 
	    BB = BB + dabs(Tol1)
          else 
	    BB = BB - dabs(Tol1)
          endif
        endif
        FB = AFunction(ftype2,a02,aux02,ct_global2,BB)
        i=i+1
      endif
	end do
    if (i >= maxIterations) error = 2
  endif
  niter = i
  BrentRoots = resultat
end Function BrentRoots

END MODULE GLOBAL
!*********************************************************************************
DOUBLE PRECISION function v0_fun(b,a)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION, INTENT(IN) :: b,a
DOUBLE PRECISION ::  slope, DCSVAL, value_left, value_right 
DOUBLE PRECISION::bvec(3),avec(3),lamb(3),value1,value2,value3,interpolate

INTEGER :: index_a,index_b,ii_b,index_bfeas,index_bfeas3,NINTV,val,nbaux,index_bfeasmax,index1,index2
EXTERNAL DCSVAL,interpolate

!uncomment to use linear interpolation
!v0_fun=interpolate(b,a,v0_matrix)
!RETURN

NINTV = nb - 1
index_a = (MAX(MIN(INT((na-1)*(a-agrid(1))/(agrid(na)-agrid(1)))+1,na-1), 1))       !closest lower grid point

!realization of a belonds to agrid
IF (a-agrid(index_a).EQ.0d0) THEN
    index1=index_bfeas4(index_a)
    val=1
    IF (((b.GE.bgrid(index1+1)).AND.(index1+1.LT.nb)).OR.(index1.EQ.0)) THEN    
        ii_b=index1+1
        nbaux =nb-ii_b+1
        v0_fun  = DCSVAL(b, nbaux-1, break_matrix(index_a,ii_b:nb), coeff_matrix(index_a,:,ii_b:nb))
        RETURN
    ELSEIF ((b.GE.bgrid(index1+1)).AND.(index1+1.EQ.nb)) THEN
        v0_fun  = v0_matrix(nb,index_a)
    ELSE
        v0_fun =-1d6 
    ENDIF
    RETURN
ENDIF

!realization of a does not belong to agrid

index1 = index_bfeas4(index_a)
index2 = index_bfeas4(index_a+1)


IF (index1.EQ.0) THEN
    if (index2.EQ.0) THEN
        ii_b=index1+1
        nbaux =nb-ii_b+1
        
        value_left  = DCSVAL(b, nbaux-1, break_matrix(index_a,ii_b:nb), coeff_matrix(index_a,:,ii_b:nb))
        value_right = DCSVAL(b, nbaux-1, break_matrix(index_a+1,ii_b:nb), coeff_matrix(index_a+1,:,ii_b:nb))

        slope = (value_right - value_left)/ (agrid(index_a+1) - agrid(index_a))
        v0_fun = value_left + slope * (a - agrid(index_a))    
        RETURN
    elseif (index2.GT.0) THEN
        if (b.GE.bgrid(index2+1)) THEN
            ii_b=index1+1
            nbaux =nb-ii_b+1
            value_left  = DCSVAL(b, nbaux-1, break_matrix(index_a,ii_b:nb), coeff_matrix(index_a,:,ii_b:nb))

            IF (index2+1.EQ.nb) THEN
                value_right = v0_matrix(nb,index_a+1)
            else                
                ii_b=index2+1
                nbaux =nb-ii_b+1        
                value_right = DCSVAL(b, nbaux-1, break_matrix(index_a+1,ii_b:nb), coeff_matrix(index_a+1,:,ii_b:nb))
            endif
            
            slope = (value_right - value_left)/ (agrid(index_a+1) - agrid(index_a))
            v0_fun = value_left + slope * (a - agrid(index_a))    
            RETURN
        elseif (b.LT.bgrid(index2+1)) THEN
            v0_fun = -1d6
            RETURN
        endif
    endif
elseif (index1.GT.0) THEN
    if (b.GE.bgrid(index1+1)) THEN
        if (index2.EQ.0) THEN
            
            IF (index1+1.EQ.nb) THEN
                value_left = v0_matrix(nb,index_a)
            else                
                ii_b=index1+1
                nbaux =nb-ii_b+1
                value_left  = DCSVAL(b, nbaux-1, break_matrix(index_a,ii_b:nb), coeff_matrix(index_a,:,ii_b:nb))
            endif           

            ii_b=index2+1
            nbaux =nb-ii_b+1        
            value_right = DCSVAL(b, nbaux-1, break_matrix(index_a+1,ii_b:nb), coeff_matrix(index_a+1,:,ii_b:nb))

            slope = (value_right - value_left)/ (agrid(index_a+1) - agrid(index_a))
            v0_fun = value_left + slope * (a - agrid(index_a))  
            RETURN
        elseif (index2.GT.0) THEN
            if (b.GE.bgrid(index2+1)) THEN
                IF (index1+1.EQ.nb) THEN
                    value_left = v0_matrix(nb,index_a)
                else                
                    ii_b=index1+1
                    nbaux =nb-ii_b+1
                    value_left  = DCSVAL(b, nbaux-1, break_matrix(index_a,ii_b:nb), coeff_matrix(index_a,:,ii_b:nb))
                endif
            
                IF (index2+1.EQ.nb) THEN
                    value_right = v0_matrix(nb,index_a+1)
                else                
                    ii_b=index2+1
                    nbaux =nb-ii_b+1        
                    value_right = DCSVAL(b, nbaux-1, break_matrix(index_a+1,ii_b:nb), coeff_matrix(index_a+1,:,ii_b:nb))
                endif
                
                slope = (value_right - value_left)/ (agrid(index_a+1) - agrid(index_a))
                v0_fun = value_left + slope * (a - agrid(index_a))    
                RETURN
            elseif (b.LT.bgrid(index2+1)) THEN
                v0_fun = -1d6
                RETURN
            endif
        endif
    elseif (b.LT.bgrid(index1+1)) THEN
        v0_fun =-1d6
        RETURN
        
    endif
endif
    
end function v0_fun
!*********************************************************************************
DOUBLE PRECISION function vaut_fun(a)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION, INTENT(IN) :: a
DOUBLE PRECISION::slope,value_left,value_right
INTEGER::index_a

index_a=(MAX(MIN(INT((na-1)*(a-agrid(1))/(agrid(na)-agrid(1)))+1,na-1), 1))

value_left  = vaut_matrix(1,index_a)
value_right = vaut_matrix(1,index_a+1)

slope = (value_right - value_left)/ (agrid(index_a+1) - agrid(index_a))

vaut_fun = value_left + slope * (a - agrid(index_a))

end function vaut_fun
!*******************************************************************************************
!FUNCTION USED TO COMPUTE THE MINIMUM INCOME FOR WHICH THE CURRENT PARTY IN POWER DOES NOT DEFAULT
!b = outstanding debt
DOUBLE PRECISION function a_fun(b,a)

USE GLOBAL

IMPLICIT NONE

REAL(8), INTENT(IN)::b,a
INTEGER :: MAXFN
DOUBLE PRECISION :: dif_fun, ERRABS, ERRREL, left, right, a_max1, a_min1old, a_min1, dif_right, dif_left
EXTERNAL dif_fun, DZBREN

DOUBLE PRECISION ::vaut_fun,v0_fun
EXTERNAL::vaut_fun,v0_fun

INTEGER::i_a,index_a,index_amin1,index_amax1,index_bfeas,index_b1,index_b2 
DOUBLE PRECISION::a1,a2,b1,b2

MAXFN=1000
ERRREL = 1D-10
ERRABS = 1D-10

a_max1 = rhoa*a + (1-rhoa)*mua + width*epsmax
a_min1 = rhoa*a + (1-rhoa)*mua + width*epsmin

dif_left  = dif_fun(a_min1,b)
dif_right = dif_fun(a_max1,b)

!1) DETERMINE WHETHER THERE IS AN INTERIOR ROOT OR NOT
!   a) IF THERE IS AN INTERIOR ROOT, SPECIFY WHETHER THE difference function IS INCREASING IN g OR NOT.
!   b) IF THERE IS NO INTERIOR ROOT, DETERMINE WHETHER THE GOV'T ALWAYS OR NEVER DEFAULTS
        
if (dif_right*dif_left<0) then  !THERE IS AN INTERIOR ROOT

    IF (dif_left.LE.-1d3+100d0) THEN
        
        if (dif_left>0) then
            a_fun = a_min1
        ELSE
            left =  a_min1
            right = a_max1
            
            CALL DZBREN (dif_fun2, ERRABS, ERRREL, left, right, MAXFN)
            a_fun = right
        endif
    else
    
        left =  a_min1
        right = a_max1
        CALL DZBREN (dif_fun2, ERRABS, ERRREL, left, right, MAXFN)
        a_fun = right
    endif
    
else if (dif_left>0) then
      a_fun = a_min1
else !dif_right<0
      a_fun = a_max1
end if

CONTAINS

DOUBLE PRECISION FUNCTION dif_fun2(x)

DOUBLE PRECISION, INTENT(IN)::x

dif_fun2 = dif_fun(x,b)

end function

end function a_fun
!!*********************************************************************************
!COMPUTE QUADRATURE POINTS AND WEIGHTS USING LEGENDRE AND GAUSSIAN QUADRATURE RULES
!THEY ARE USED TO COMPUTE NUMERICAL INTEGRALS
subroutine quadrature

USE GLOBAL

IMPLICIT NONE

INTEGER :: N_v, N_q, N, IWEIGH, IWEIGH4, NFIX
parameter(N_v = nquad_v, N_q = nquad_q, N = neps)
DOUBLE PRECISION:: ALPHA1, BETA1, QW_v(1:N_v), QX_v(1:N_v), QW_q(1:N_q), QX_q(1:N_q),QW(1:N), QX(1:N), QXFIX(2)
PARAMETER(ALPHA1=0, BETA1=0, IWEIGH=1, IWEIGH4=4)
external DGQRUL

!quad_w = weights using Gauss legendre quadrature rule
!quad_x = points using Gauss legendre quadrature rule
!quad_w_hermite = weights using Gauss Hermite quadrature rule
!quad_x_hermite = points using Gauss Hermite quadrature rule

NFIX = 0
CALL DGQRUL(N_v,IWEIGH, ALPHA1, BETA1, NFIX, QXFIX, QX_v, QW_v)

quad_w_v=QW_v
quad_x_v=QX_v

CALL DGQRUL(N_q,IWEIGH, ALPHA1, BETA1, NFIX, QXFIX, QX_q, QW_q)

quad_w_q=QW_q
quad_x_q=QX_q

CALL DGQRUL(N,IWEIGH4, ALPHA1, BETA1, NFIX, QXFIX, QX, QW)
quad_x_hermite = QX
quad_w_hermite = QW

end subroutine
!*********************************************************************************
DOUBLE PRECISION function interpolate(b,a,matrix)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION, INTENT(IN) :: b,a, matrix(nb,na)
DOUBLE PRECISION ::  slope, weight, acum, ratio_b,ratio_a
INTEGER :: index_b,index_a,i_b,i_a

index_a = (MAX(MIN(INT((na-1)*(a-agrid(1))/(agrid(na)-agrid(1)))+1,na-1), 1))
index_b = (MAX(MIN(INT((nb-1)*(b-bgrid(1))/(bgrid(nb)-bgrid(1)))+1,nb-1), 1))

ratio_b = (b - bgrid(index_b)) / (bgrid(index_b+1) - bgrid(index_b))
ratio_a = (a - agrid(index_a)) / (agrid(index_a+1) - agrid(index_a))
acum=0

do i_b=0,1
   do i_a=0,1
        weight = ((1-i_b)*(1 - ratio_b) + i_b*ratio_b )* ((1-i_a)*(1 - ratio_a) + i_a*ratio_a )
        acum = acum + matrix(index_b + i_b, index_a + i_a)*weight
   end do
end do

interpolate = acum
end
!*********************************************************************************
DOUBLE PRECISION function interpolate3(a,b,matrix)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION, INTENT(IN) :: b,a, matrix(na,nb)
DOUBLE PRECISION ::  slope, weight, acum, ratio_b,ratio_a
INTEGER :: index_b,index_a,i_b,i_a

index_a = (MAX(MIN(INT((na-1)*(a-agrid(1))/(agrid(na)-agrid(1)))+1,na-1), 1))
index_b = (MAX(MIN(INT((nb-1)*(b-bgrid(1))/(bgrid(nb)-bgrid(1)))+1,nb-1), 1))

ratio_b = (b - bgrid(index_b)) / (bgrid(index_b+1) - bgrid(index_b))
ratio_a = (a - agrid(index_a)) / (agrid(index_a+1) - agrid(index_a))
acum=0

do i_b=0,1
   do i_a=0,1
        weight = ((1-i_b)*(1 - ratio_b) + i_b*ratio_b )* ((1-i_a)*(1 - ratio_a) + i_a*ratio_a )
        acum = acum + matrix(index_a + i_a, index_b + i_b)*weight
   end do
end do

interpolate3 = acum
end
!*********************************************************************************
!SPECIFY GRID VALUES
subroutine compute_bgrid

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION :: dos
INTEGER :: i0(1),i

do i=1,nb
   bgrid(i) = bmin + (bmax - bmin) * REAL((i-1d0) / (nb - 1d0))
end do

i0 = LOCATE(bgrid,0d0)
i_bzero = i0(1)
bgrid(i_bzero)=0d0


OPEN(1,FILE='graphs//bgrid.txt')

DO i = 1,nb
    WRITE(1,*) bgrid(i)
ENDDO
    
CLOSE(1)

end subroutine
!*********************************************************************************
DOUBLE PRECISION function dif_fun(a,b)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION,INTENT(IN)::a,b
DOUBLE PRECISION::v0_fun,vaut_fun,dif_fun_temp
EXTERNAL v0_fun,vaut_fun

dif_fun = v0_fun(b,a)-vaut_fun(a)
dif_fun_temp = dif_fun

end function dif_fun
!!*********************************************************************************
REAL(8) FUNCTION udef(a)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::a

IF (choice_defcost.EQ.1) THEN
    udef = 0d0
ELSEIF (choice_defcost.EQ.2) THEN
    udef = MAX(0d0,psi1+psi2*DLOG(DEXP(a)))
ENDIF

END FUNCTION udef
!*********************************************************************************
!FUNCTION THAT COMPUTES THE PRICES ASKED FOR A COUNTRY BOND
!b_next = future debt
!y = current outcome
DOUBLE PRECISION function q_fun(b_next,a)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION, INTENT(IN)::b_next,a
DOUBLE PRECISION :: prob_no_def_tomorrow, a_current_type, DNORDF, a_fun, standarized_inf, standarized_sup,&
                    a_threshold, a_min1, a_max1, exp_q, scalar, a_next, q_paid_fun, acum_int, &
                    cdf, cdf_threshold, DNORIN
INTEGER :: i, NINTV
EXTERNAL DNORDF, a_fun, q_paid_fun, DNORIN

!a_current_type = STANDARIZED INCOME THRESHOLDS IF THE TYPE DOES NOT CHANGE TOMORROW

if (psi1.GE.500d0) then
   q_fun = payoff/ (lambda + r)
    RETURN
ENDIF

if (b_next >=0) then
   q_fun = payoff/ (lambda + r)
else

   a_threshold = a_fun(b_next,a)
   a_current_type = (a_threshold - rhoa * a - (1-rhoa)*mua) / sigmaa

   standarized_inf = epsgrid(1) / sigmaa
   standarized_sup = epsgrid(neps) / sigmaa

   a_min1 = (1-rhoa)*mua + rhoa* a - 4*sigmaa 
   a_max1 = (1-rhoa)*mua + rhoa* a + 4*sigmaa 

   scalar = 1d+0 / SQRT(pi_number)
   NINTV = ncdf - 1
   exp_q = 0d+0
   acum_int = 0d+0
      
      if (DNORDF(a_current_type) <= DNORDF(standarized_inf)+1d-10) then !threshold tomorrow is low
         prob_no_def_tomorrow = 1d+0
         do i=1,nquad_q
             cdf = 0.5*(quad_x_q(i) + 1)*(cdfmax - cdfmin) + cdfmin
             a_next  = (1-rhoa)*mua + rhoa* a + DNORIN(cdf) * sigmaa
             exp_q = exp_q + (1d+0 - lambda) * q_paid_fun(b_next,a_next) * 0.5* (cdfmax - cdfmin)*quad_w_q(i)
         end do
         exp_q = exp_q / (cdfmax-cdfmin)  
   elseif (DNORDF(a_current_type) >= DNORDF(standarized_sup)-1d-10) then !threshold tomorrow is high
         prob_no_def_tomorrow = 0d+0
         exp_q = 0d+0
   else
         cdf_threshold = DNORDF(a_current_type)
         prob_no_def_tomorrow = (1d+0 - cdf_threshold)
         do i=1,nquad_q
             !compute integral for y > y_threshold
             cdf = 0.5*(quad_x_q(i) + 1)*(cdfmax- cdf_threshold) + cdf_threshold
             a_next  = (1-rhoa)*mua + rhoa* a + DNORIN(cdf) * sigmaa
             exp_q = exp_q + (1d+0 - lambda) * q_paid_fun(b_next,a_next) * 0.5* (cdfmax - cdf_threshold)*quad_w_q(i)
         end do
         exp_q = exp_q * prob_no_def_tomorrow / (cdfmax - cdf_threshold)  
   end if

q_fun = MAX(prob_no_def_tomorrow*payoff / (1d+0 + r) + exp_q/(1d+0 + r), 0d+0)

end if

end
!*********************************************************************************
DOUBLE PRECISION function q_paid_fun(b,a)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION, INTENT(IN)::b,a
DOUBLE PRECISION ::  slope,DCSVAL,value_left, value_right
INTEGER :: index_a,NINTV,i_b
EXTERNAL DCSVAL

if (psi1.GE.500d0) then
   q_paid_fun = payoff/ (lambda + r)
   RETURN
ENDIF


NINTV = nb - 1
index_a = (MAX(MIN(INT((na-1)*(a-agrid(1))/(agrid(na)-agrid(1)))+1,na-1), 1))
    
value_left  = DCSVAL(b, NINTV, break_matrix_q(index_a,:), coeff_matrix_q(index_a,:,:))
value_right = DCSVAL(b, NINTV, break_matrix_q(index_a+1,:), coeff_matrix_q(index_a+1,:,:))

slope = (value_right - value_left)/ (agrid(index_a+1) - agrid(index_a))

q_paid_fun = value_left + slope * (a- agrid(index_a))

end
!***************************************************************************
DOUBLE PRECISION function q_fun_base(b,a)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION, INTENT(IN)::b,a
DOUBLE PRECISION ::  slope,DCSVAL,value_left, value_right
INTEGER :: index_a,NINTV,index_bfeas,i_b
REAL(8)::interpolate2,q_fun
EXTERNAL::q_fun

if (psi1.GE.500d0) then
   q_fun_base = payoff/ (lambda + r)
   RETURN
ENDIF

q_fun_base = interpolate2(b,a,q1base)

q_fun_base=MIN(q_fun_base,payoff/ (lambda + r))
q_fun_base=MAX(q_fun_base,0d0)

RETURN
end
!**********************************************************************
DOUBLE PRECISION function interpolate2(b,a,matrix)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION, INTENT(IN) :: b,a, matrix(na,501)
DOUBLE PRECISION ::  slope, weight, acum, ratio_b,ratio_a
INTEGER :: index_b,index_a,i_b,i_a,nb3

nb3=501

index_a = (MAX(MIN(INT((na-1)*(a-agrid(1))/(agrid(na)-agrid(1)))+1,na-1), 1))
index_b = (MAX(MIN(INT((nb3-1)*(b-bgridbase(1))/(bgridbase(nb3)-bgridbase(1)))+1,nb3-1), 1))

ratio_b = (b - bgridbase(index_b)) / (bgridbase(index_b+1) - bgridbase(index_b))
ratio_a = (a - agrid(index_a)) / (agrid(index_a+1) - agrid(index_a))
acum=0

do i_b=0,1
   do i_a=0,1
        weight = ((1-i_b)*(1 - ratio_b) + i_b*ratio_b )* ((1-i_a)*(1 - ratio_a) + i_a*ratio_a )
        acum = acum + matrix(index_a + i_a, index_b + i_b )*weight
   end do
end do

interpolate2 = acum
end
!*********************************************************************************